#ifndef COLLISION_H
#define COLLISION_H

#include <vector>
#include <string>
#include <fstream>


#include "../Tool/Macros.h"
#include "../Input/PicParams.h"
#include "../Species/Species.h"
#include "../Diagnostic/Diagnostic.h"
#include "../Diagnostic/Diagnostic1D.h"
#include "../Diagnostic/Diagnostic2D.h"
#include "../Tool/Log.h"

class Diagnostic;

using namespace std;

class Collision
{

public:

    Collision(PicParams* params_, CollisionParameter collision_param_, int collision_number_):
        params(params_),
        collision_param(collision_param_),
        collision_number(collision_number_)
    {

    };

    virtual ~Collision(){};

    //non-relativistic case
    virtual void collide(PhysicalField* fields, vector<Species*>& vecSpecies, Diagnostic* diag, int itime){};

    virtual void read_cross_section()
    {
        double energy, cross_section;
        crossSection.resize(2);

        //inFile.open(collision_param.cs_file_name.c_str());
        ifstream inFile(collision_param.cs_file_name.c_str());
        if(!inFile)
        {
            log_error<<"no cross section file: "<<collision_param.cs_file_name;
        }

        while(inFile >> energy && inFile >> cross_section)
        {
            crossSection[0].push_back(energy);
            crossSection[1].push_back(cross_section);
            //crossSection[1].push_back(cross_section);
        }
        inFile.close();
    };

    //interplate cross section with energy (eV)
    inline double interplate_cross_section(double energy)
    {
        int low = 0;
        int high = crossSection[0].size() - 1;
        int mid = 0;

        if(energy <= crossSection[0][low] )
        {
            return 0.0;
        }
        else if(energy >= crossSection[0][high])
        {
            low = high -1;
            double dEnergy_inv = 1.0 / (crossSection[0][low] - crossSection[0][high]);
            return crossSection[1][high] * (crossSection[0][low] - energy) * dEnergy_inv +
                   crossSection[1][low] * (energy - crossSection[0][high]) * dEnergy_inv;
        }

        while(low <= high){
            mid = (low + high) / 2;
            if(crossSection[0][mid] < energy){
                low = mid + 1;
            }
            else if(crossSection[0][mid] > energy){
                high = mid - 1;
            }
            else {
                return crossSection[1][mid];
            }
        }
        //now low is 1 larger than high
        double dEnergy_inv = 1.0 / (crossSection[0][low] - crossSection[0][high]);
        return crossSection[1][high] * (crossSection[0][low] - energy) * dEnergy_inv +
               crossSection[1][low] * (energy - crossSection[0][high]) * dEnergy_inv;
    };

    //Simple method: Calculate electron scattered velocity for electron-neutral collision
    //the method is eqution (11) from the ref: a Monte Carlo collision model for the particle in cell method: applications to
    //argon and oxygen discharges.
    //and the code is transformed from C.F. Sang's fortran code
    inline void calculate_scatter_velocity_simple(double v_magnitude, double mass1, double mass2,
                                    vector<double>& momentum_unit, vector<double>& momentum_temp)
    {
        double up1, up2, up3;
        double r11, r12, r13, r21, r22, r23, r31, r32, r33;
        double mag;

        double ra = (double)rand() / RAND_MAX;
        double costheta = 1.0 - 2.0 * ra;
        double sintheta = sqrt(1.0 - abs(costheta * costheta) );

        ra = (double)rand() / RAND_MAX;
        double pi = 3.1415926;
        double phi = 2.0 * pi * ra;
        double cosphi = cos(phi);
        double sinphi = sin(phi);

        double ve=v_magnitude*sqrt(1.0-2.0*mass1*(1.0-costheta)/mass2);

        r13 = momentum_unit[0];
        r23 = momentum_unit[1];
        r33 = momentum_unit[2];
        if(r33 == 1.0 ){
            up1= 0.;
            up2= 1.;
            up3= 0.;
        }
        else{
            up1= 0.;
            up2= 0.;
            up3= 1.;
        }

        r12 = r23 * up3 - r33 * up2;
        r22 = r33 * up1 - r13 * up3;
        r32 = r13 * up2 - r23 * up1;
        mag = sqrt(r12 * r12 + r22 * r22 + r32 * r32);
        r12 = r12 / mag;
        r22 = r22 / mag;
        r32 = r32 / mag;
        r11 = r22 * r33 - r32 * r23;
        r21 = r32 * r13 - r12 * r33;
        r31 = r12 * r23 - r22 * r13;
        momentum_temp[0] = ve * (r11 * sintheta * cosphi + r12 * sintheta * sinphi + r13 * costheta);
        momentum_temp[1] = ve * (r21 * sintheta * cosphi + r22 * sintheta * sinphi + r23 * costheta);
        momentum_temp[2] = ve * (r31 * sintheta * cosphi + r32 * sintheta * sinphi + r33 * costheta);
    };


    //Complex method: Calculate electron scattered velocity for electron-neutral collision
    inline void calculate_scatter_velocity(double ke, double v_magnitude, double mass1, double mass2,
                                        vector<double>& momentum_unit, vector<double>& momentum_temp)
    {
        double up1, up2, up3;
        double r11, r12, r13, r21, r22, r23, r31, r32, r33;
        double mag;

        double ra = (double)rand() / RAND_MAX;
        double cosX = ( 2.0 + ke - 2.0 * pow(1.0+ke, ra) ) / ke;
        //double cosX = 1.0 - 2.0 * ra;
        double sinX = sqrt(1.0 - abs(cosX * cosX) );
        while( cosX*cosX > 1.0 || std::isnan(sinX) )
        {
            ra = (double)rand() / RAND_MAX;
            cosX = ( 2.0 + ke - 2.0 * pow(1.0+ke, ra) ) / ke;
            //cosX = 1.0 - 2.0 * ra;
            sinX = sqrt(1.0 - abs(cosX * cosX) );
        }

        ra = (double)rand() / RAND_MAX;
        double pi = 3.1415926;
        double phi = 2.0 * pi * ra;
        double cosphi = cos(phi);
        double sinphi = sin(phi);
        double sinTheta_inv;

        double ve=v_magnitude*sqrt(1.0-2.0*mass1*(1.0-cosX)/mass2);

        r11 = momentum_unit[0];
        r12 = momentum_unit[1];
        r13 = momentum_unit[2];

        if(r11 == 1.0 || r11 == -1.0)
        {
            up1= 0.;
            up2= 1.;
            up3= 0.;
        }
        else
        {
            up1= 1.;
            up2= 0.;
            up3= 0.;
        }

        //double cosTheta = r11 * up1 + r12 * up2 + r13 * up3;
        //sinTheta_inv = 1.0 / sqrt(1.0 - cosTheta*cosTheta);

        r21 = ( r12 * up3 - r13 * up2 );
        r22 = ( r13 * up1 - r11 * up3 );
        r23 = ( r11 * up2 - r12 * up1 );

        //the code below is not consistent with the ref, but this is right.
        //and the ref may be wrong in somewhere???
        mag = sqrt( r21*r21 + r22*r22 + r23*r23 );
        r21 /= mag;
        r22 /= mag;
        r23 /= mag;

        r31 = ( r22 * r13 - r23 * r12 );
        r32 = ( r23 * r11 - r21 * r13 );
        r33 = ( r21 * r12 - r22 * r11 );


        momentum_temp[0] = ve * (r11 * cosX + r21 * sinX * sinphi + r31 * sinX*cosphi);
        momentum_temp[1] = ve * (r12 * cosX + r22 * sinX * sinphi + r32 * sinX*cosphi);
        momentum_temp[2] = ve * (r13 * cosX + r23 * sinX * sinphi + r33 * sinX*cosphi);

    };


    vector< vector<double> > crossSection;

    PicParams* params;

    const CollisionParameter collision_param;

    const int collision_number;

private:



};


#endif
